|
\[\textbf{Lineu Alberto Cavazani de Freitas}\] \[\textbf{Prof. Cesar Augusto Taconeli}\] \[\textbf{Prof. José Luiz Padilha da Silva}\]
Análise Comportamental de Ovelhas Submetidas a Intervenção Humana usando GAMLSS
Material para Paper: Número de Mudanças de Postura de Orelha e Proporção com as Orelhas em Posição Neutra
1 Experimento
O conjunto de dados é proveniente de um experimento sobre o comportamento de ovelhas, conduzido na fazenda experimental INRA La Fage, Roqueford, França, em setembro de 2015 com o objetivo de verificar o efeito de linhagem genética, escovação e isolamento nas respostas comportamentais dos animais (TAMIOSO et al., 2017). Na ocasião, vinte ovelhas classificadas como reativas ou não reativas ao isolamento social temporário foram submetidas à escovação por um humano familiar. As ovelhas tinham 15 meses de idade, eram não gestantes e não amamentavam quando foram observadas.
|
O experimento foi conduzido em três sessões experimentais: na primeira tinha-se uma grade de metal separando o animal testado dos demais animais, sem distância entre eles. Na segunda havia duas grades de metal separando os animais a uma distância de 1,7 metros, ou seja, foi imposta a condição de isolamento social. E na terceira sessão, os animais voltaram a ser separados por apenas uma grade:
|
As sessões de testes ocorreram dois dias após a fase de adaptação dos animais ao equipamento e aos humanos e, em cada sessão, as ovelhas foram observadas em 3 momentos distintos: fase de pré escovação, com duração de 2 minutos e 30 segundos; fase de escovação, com duração de 3 minutos; e pós escovação, com duração de 2 minutos e 30 segundos. O delineamento para um animal pode ser represntado da seguinte forma:
|
Os dados coletados dizem respeito ao número de mudanças de postura dos animais e à proporção do tempo em que os animais permaneceram em determinadas posturas, tratando-se então, de um conjunto de dados com múltiplas respostas em que não há observações independentes, já que cada animal contribui com nove medidas. Portanto, há a necessidade de incorporar as correlações entre as medidas num mesmo animal e do animal dentro de cada sessão experimental, além da correlação entre as respostas.
Foram avaliados os efeitos de:
Sessão: Fator de 3 níveis que indica a sessão experimental (Se1, Se2, Se3).
Momento: Fator de 3 níveis em que indica o momento experimental (antes, durante, depois).
Linhagem: Fator de 2 níveis que classifica os animais como reativos ou não reativos ao isolamento social temporário.
2 Análise Exploratória
2.1 Dados
setwd("~/Área de Trabalho/PESQUISA/IC_nova/1. Dados")
######################################################################
## Pacotes
library(ggplot2)
library(gridExtra)
library(gamlss)
######################################################################
## Dados
### Contagem
count <- read.csv2("DataCont.csv", header = T, sep = ',')
## Alterando os contrastes default para os contrastes de interesse
count$sessaoR <- count$sessao
count$sessaoR <- relevel(count$sessaoR, ref = 'Se3')
contrasts(count$sessaoR) <- 'contr.helmert'
contrasts(count$sessaoR)[,1] <- contrasts(count$sessaoR)[,1]/2
contrasts(count$sessaoR)[,2] <- contrasts(count$sessaoR)[,2]/3
## Criando um fator combinando animal e sessao
count$animals <- factor(count$animal):count$sessao
### Proporção
prop <- read.csv2("DataProp.csv", header = T, sep = ';')
prop <- prop[,-c(8)]
## Alterando os contrastes default para os contrastes de interesse
prop$sessaoR <- prop$sessao
prop$sessaoR <- relevel(prop$sessaoR, ref = 'Se3')
contrasts(prop$sessaoR) <- 'contr.helmert'
contrasts(prop$sessaoR)[,1] <- contrasts(prop$sessaoR)[,1]/2
contrasts(prop$sessaoR)[,2] <- contrasts(prop$sessaoR)[,2]/3
## Criando um fator combinando animal e sessao
prop$animals <- factor(prop$animal):prop$sessao
## Transformação da resposta para que seja possível utilizar
## a distribuição beta
prop$neutra <- 1-prop$resporelha2
#sum(prop$neutra == 1)
#sum(prop$neutra == 0)
prop$neutra.be <- ((prop$neutra*(180-1))+(0.5))/180
#sum(prop$neutra.be == 1)
#sum(prop$neutra.be == 0)2.2 Histograma e Box-plots
2.3 Gráficos de Perfil
3 Modelos para o Número de Mudanças de Postura de Orelha
Para o ajuste dos modelos de regressão foram consideradas as distribuições:
- PO: Poisson
- NBI:Binomial Negativa Tipo I
- ZIP: Poisson Inflacionada de Zeros
- ZINBI: Binomial Negativa Inflacionada de Zeros
- ZAP: Poisson Zero Ajustada
- ZANBI: Binomial Negativa Zero Ajustada
Para cada uma das 7 distribuições ajustou-se o modelo no qual era considerado sessão, momento, linhagem, as interações duas a duas e efeitos aleatórios a nível de animal e animal dentro de sessão para a média.
Nas distribuições com parâmetros referentes à inflação em 0 também foram acrescentadas as mesmas covariáveis, com exceção dos efeitos aleatórios.
Todas análises foram feitas utilizando o pacote gamlss do software estatístico R.
modelo_norelha <- function(familia){
gamlss(norelha ~ sessaoR + tempo + linhagem
+ (sessaoR + tempo + linhagem)^2 +
re(random = list(animal=~1, animals=~1), method = 'REML'),
nu.formula = norelha ~ sessaoR + tempo +
linhagem+ (sessaoR + tempo + linhagem)^2,
data = count, method = RS(50),
family = familia)}
mPO <- modelo_norelha(familia = PO())
mNBI <- modelo_norelha(familia = NBI())
mZIP <- modelo_norelha(familia = ZIP())
mZINBI <- modelo_norelha(familia = ZINBI())
mZAP <- modelo_norelha(familia = ZAP())
mZANBI <- modelo_norelha(familia = ZANBI())3.1 Seleção da Família Via Medidas de Qualidade
A escolha da distribuição utilizada baseou-se nos valores observados para as medidas AIC, BIC, Deviance, Verossimilhança, graus de liberdade e por fim, verificou-se se houve convergência.
## Modelo AIC BIC Deviance logLik df Conv Iter
## 1 PO 1191.028 1344.848 1094.6787 -547.3393 48.17479 TRUE 5
## 2 NBI 1054.742 1206.695 959.5615 -479.7808 47.59021 TRUE 7
## 3 ZIP 1127.398 1315.057 1009.8519 -504.9260 58.77281 TRUE 17
## 4 ZINBI 2488.164 2687.903 2363.0517 -1181.5259 62.55619 FALSE 50
## 5 ZAP 1137.740 1324.007 1021.0663 -510.5332 58.33675 TRUE 12
## 6 ZANBI 1046.170 1247.390 920.1308 -460.0654 63.01976 TRUE 10A distribuição selecionada foi a Binomial Negativa Zero Ajustada.
3.2 Reajustes
Considerando o modelo ZANBI, foram realizados 3 reajustes para o modelo original:
- Retirando interações de \(\nu\)
- Retirando interações de \(\mu\) e \(\nu\)
- Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)
Todos os modelos reajustados são modelos encaixados ao modelo original, sendo assim utilizou-se o teste de razão de verossimilhanças para selecionar o menor modelo possível que não difere estatísticamente do original.
3.2.1 Retirando interações de \(\nu\)
mZANBIr1 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
(sessaoR + tempo + linhagem)^2 +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
nu.formula = norelha ~ sessaoR + tempo + linhagem,
family = ZANBI(), data = count)3.2.2 Retirando interações de \(\mu\) e \(\nu\)
mZANBIr2 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
nu.formula = norelha ~ sessaoR + tempo + linhagem,
family = ZANBI(), data = count)3.2.3 Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)
mZANBIr3 <- gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
family = ZANBI(), data = count)3.2.4 Medidas de qualidade e TRV
## Modelo AIC BIC Deviance logLik df Conv P-val
## 1 Completo 1046.170 1247.390 920.1308 -460.0654 63.01976 TRUE -
## 2 Reajuste 1 1036.359 1212.034 926.3192 -463.1596 55.01976 TRUE 0.63
## 3 Reajuste 2 1022.671 1174.258 927.7194 -463.8597 47.47563 TRUE 0.95
## 4 Reajuste 3 1037.248 1172.871 952.2970 -476.1485 42.47563 TRUE 0.05Conclusões
O modelo sem interações para nu não difere do modelo completo.
O modelo sem interações para mu e nu não difere do modelo completo.
O modelo sem interações para mu e nenhuma covariável para nu difere do modelo com covariáveis para nu.
Sendo assim, as covariáveis para mu e nu devem ser mantidas, porém sem as interações.
4 Modelos para a Proporção do Tempo com as Orelhas em Posição Neutra
Para o ajuste dos modelos de regressão foram consideradas as distribuições:
- BE: Beta
- BEINF:Beta Inflacionada
Para cada uma das distribuições ajustou-se o modelo no qual era considerado sessão, momento, linhagem, as interações duas a duas e efeitos aleatórios a nível de animal e animal dentro de sessão para a média.
Na distribuição Beta Inflacionada acrescentou-se as mesmas covariáveis, com exceção dos efeitos aleatórios no parâmetro que controla a inflação em 0.
m1.be <- gamlss(neutra.be ~ sessaoR + tempo + linhagem
+ (sessaoR + tempo + linhagem)^2 +
re(random = list(animal=~1, animals=~1), method = 'REML'),
data = prop, family = BE())
m1.beinf <- gamlss(neutra ~ sessaoR + tempo + linhagem
+ (sessaoR + tempo + linhagem)^2 +
re(random = list(animal=~1, animals=~1), method = 'REML'),
nu.formula = ~ sessaoR + tempo + linhagem,
data = prop, family = BEINF())4.1 Seleção da Família Via Medidas de Qualidade
A escolha da distribuição utilizada baseou-se na análise gráfica dos resíduos.
A análise gráfica evidencia um comportamento satisfatório dos resíduos no modelo com distribuição Beta Inflacionada.
4.2 Reajustes
Considerando o modelo BEINF, foram realizados 2 reajustes para o modelo original:
- Retirando interações de \(\mu\)
- Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)
Os modelos reajustados são modelos encaixados ao modelo original, sendo assim utilizou-se o teste de razão de verossimilhanças para selecionar o menor modelo possível que não difere estatísticamente do original.
4.2.1 Retirando interações de \(\mu\)
m1.beinf.r1 <- gamlss(formula = neutra ~ sessaoR + tempo + linhagem +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
nu.formula = ~ sessaoR + tempo + linhagem,
family = BEINF(), data = prop)4.2.2 Retirando interações de \(\mu\) e todas as covariáveis de \(\nu\)
m1.beinf.r2 <- gamlss(formula = neutra ~ sessaoR + tempo + linhagem +
re(random = list(animal = ~1, animals = ~1), method = 'REML'),
family = BEINF(), data = prop)4.2.3 Medidas de qualidade e TRV
## Modelo AIC BIC Deviance logLik df P-val
## 1 Completo 98.71873 286.7912 -19.085821 9.542911 58.90228 -
## 2 Reajuste 1 98.34053 249.2352 3.823311 -1.911655 47.25861 0.02
## 3 Reajuste 2 139.74139 274.6713 55.224177 -27.612088 42.25861 05 Análise de Resíduos
Parte importante da etapa de ajuste de modelos de regressão é a análise de resíduos. A análise de resíduos em modelos da classe GAMLSS é baseada em resíduos quantílicos aleatorizados, estes resíduos são de interpretação mais simples visto que, sob um modelo bem ajustado, possuem distribuição Normal.
6 Resumo dos modelos finais
6.1 Número de Mudanças de Postura de Orelha
## ******************************************************************
## Family: c("ZANBI", "Zero Altered Negative binomial type I")
##
## Call: gamlss(formula = norelha ~ sessaoR + tempo + linhagem +
## re(random = list(animal = ~1, animals = ~1), method = "REML"),
## nu.formula = norelha ~ sessaoR + tempo + linhagem,
## family = ZANBI(), data = count)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.49220 0.07710 32.324 < 2e-16 ***
## sessaoR1 0.42218 0.10556 3.999 0.000105 ***
## sessaoR2 -0.16399 0.09251 -1.773 0.078580 .
## tempoDepois -0.43246 0.09736 -4.442 1.86e-05 ***
## tempoDurante -1.05623 0.11387 -9.276 4.36e-16 ***
## linhagems+ 0.28393 0.08580 3.309 0.001205 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9154 0.2149 -8.912 3.43e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: logit
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.1780 1.0714 -3.899 0.000153 ***
## sessaoR1 -2.6613 1.0911 -2.439 0.016049 *
## sessaoR2 0.4330 0.7312 0.592 0.554710
## tempoDepois 1.4892 1.1487 1.296 0.197069
## tempoDurante 2.7598 1.0831 2.548 0.011973 *
## linhagems+ -1.1191 0.6379 -1.754 0.081693 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 180
## Degrees of Freedom for the fit: 47.47563
## Residual Deg. of Freedom: 132.5244
## at cycle: 14
##
## Global Deviance: 927.7194
## AIC: 1022.671
## SBC: 1174.258
## ******************************************************************6.2 Proporção do Tempo com as Orelhas em Posição Neutra
## ******************************************************************
## Family: c("BEINF", "Beta Inflated")
##
## Call: gamlss(formula = neutra ~ sessaoR + tempo + linhagem +
## (sessaoR + tempo + linhagem)^2 + re(random = list(animal = ~1,
## animals = ~1), method = "REML"), nu.formula = ~sessaoR +
## tempo + linhagem, family = BEINF(), data = prop)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: logit
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.83449 0.25094 -7.310 3.13e-11 ***
## sessaoR1 -0.01893 0.49691 -0.038 0.96968
## sessaoR2 -0.92998 0.43303 -2.148 0.03374 *
## tempoDepois 0.12342 0.32152 0.384 0.70175
## tempoDurante 0.39362 0.34255 1.149 0.25278
## linhagems+ 0.28003 0.30890 0.907 0.36645
## sessaoR1:tempoDepois 0.75939 0.62393 1.217 0.22593
## sessaoR2:tempoDepois 0.20267 0.47892 0.423 0.67292
## sessaoR1:tempoDurante 0.72923 0.56750 1.285 0.20125
## sessaoR2:tempoDurante -1.29466 0.53008 -2.442 0.01604 *
## sessaoR1:linhagems+ 1.17052 0.50606 2.313 0.02241 *
## sessaoR2:linhagems+ 0.62529 0.43144 1.449 0.14984
## tempoDepois:linhagems+ 0.51738 0.42965 1.204 0.23086
## tempoDurante:linhagems+ 1.40139 0.43795 3.200 0.00176 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: logit
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.36897 0.09445 -3.906 0.000155 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: log
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6365 0.3515 -1.811 0.0726 .
## sessaoR1 -2.8428 0.5109 -5.564 1.61e-07 ***
## sessaoR2 0.4536 0.3689 1.230 0.2212
## tempoDepois 0.3786 0.4366 0.867 0.3876
## tempoDurante 1.0164 0.4408 2.306 0.0228 *
## linhagems+ -0.9210 0.3627 -2.540 0.0124 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.673 1.005 -4.651 8.49e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 180
## Degrees of Freedom for the fit: 58.90228
## Residual Deg. of Freedom: 121.0977
## at cycle: 12
##
## Global Deviance: -19.08582
## AIC: 98.71873
## SBC: 286.7912
## ******************************************************************TAMIOSO, P. R.; BOISSY, A.; BOIVIN, X.; CHANDEZE, H.; ANDANSON, S.; DELVAL, E.; HAZARD, D.; TACONELI, C. A.; MOLENTO, C. F. M. Does emotional reactivity influence behavioral and cardiac responses of ewes submitted to brushing? Behavioural Processes, p. np, 2017.
|
|
|
|
|